Gaussian phase-space representations for fermions 
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We introduce a positive phase-space representation for fermions, using the most general possi- 
ble multi-mode Gaussian operator basis. The representation generalizes previous bosonic quantum 
phase-space methods to Fermi systems. We derive equivalences between quantum and stochastic 
moments, as well as operator correspondences that map quantum operator evolution onto stochastic 
processes in phase space. The representation thus enables first-principles quantum dynamical or 
equilibrium calculations in many-body Fermi systems. Potential applications are to strongly inter- 
acting and correlated Fermi gases, including coherent behaviour in open systems and nanostructures 
described by master equations. Examples of an ideal gas and the Hubbard model are given, as well 
as a generic open system, in order to illustrate these ideas. 
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I. INTRODUCTION 

The study of strongly correlated Fermi gases is one of 
the most active areas in modern condensed matter and 
AMO physics. In quantum degenerate electron gases, im- 
provements in condensed matter materials have led to so- 
phisticated experiments, typically in reduced dimensional 
environments. Many interesting quantum phenomena are 
observed in these systems, including such features as the 
quantum Hall effect, metal-insulator phase-transitions 2 , 
high T-c superconductors^, and single electron gates in 
nanostructures^. 

Recently, pioneering experiments in strongly- 
interacting ultra-cold Fermi gases have opened up 
novel experiments of unprecedented simplicity and 
precision, both in the BEC-BCS cross-over regime^, and 
in lattices 6 . The underlying atomic interactions are ex- 
tremely well-understood, and the dynamics, interactions 
and geometry are all highly adaptable. Measurement 
techniques are also rapidly improving, with direct 
measurements of collective modest, thermodynamic 
properties^, vortices^, and even momentum correlations 
being recently reported^. 

This situation provides a substantial opportunity to 
develop and test first-principles theoretical methods for 
the investigation of correlations and dynamical effects in 
quantum degenerate Fermi gases. To this end, we intro- 
duce a generalised phase-space representation for corre- 
lated fermionic systems. The representation is based on a 
Gaussian operator basis for fermionic density operators. 
Like the analogous basis for bosonsii, the fermionic oper- 
ator basis enables the representation of arbitrary physical 
density operators as a positive distribution over a phase 
space. This representation allows quantum evolution, ei- 
ther in real time or in inverse temperature, to be viewed 
as a stochastic evolution of covariances or Green's func- 
tions. 

Phase-space methods based on coherent states^ have 
long been used for bosonic systems, with great success. 
These approaches include the Wigner function^, the Q- 
functionia, as well as the well-known Glauber-Sudarshan 
P- function^, and its generalizationsi^ii 7 .. The early 



methods based on classical phase spaces were later gen- 
eralized to give the positive-P distribution^, which has 
proved a successful way to simulate quantum many-body 
systems from first principles 19 . This method reduces 
quantum dynamics to the time evolution of a positive 
distribution on an over-complete basis set of coherent 
state projection operators, which are special cases of the 
bosonic Gaussian operators. Applications have been to 
quantum statistics of lasers 2 ^ superfluorescenceSI, para- 
metric amplifiersiSi 2 ^, quantum solitons 2 ^, as well as 
quantum dynamics 2 ^ and thermal correlations 2 ^ in Bose- 
Einstein condensates. 

Fermionic phase-space representations are relevant to 
a long-standing problem in theoretical physics, which is 
the sign problem that occurs in many-body fermionic 
wave-functionsS&SL 2 ^. There are many different approx- 
imate techniques that can be used, but the intention of 
this paper is to establish fundamentally exact procedures 
to treat the Fermi sign problem. As shown inS, the 
Gaussian method can be applied to the difficult case of 
the repulsive Hubbard model^S. Here we concentrate on 
the foundational issues of the Gaussian representation 
method, presenting the general identities required to ap- 
ply the method to a wide range of problems in fermionic 
many-body physics, including both ultra-cold atomic and 
condensed matter systems. 

To proceed, we make use of three important results, 
proved elsewhere^!: 

• the Gaussian fermion operators form a complete 
basis for any physical density operator, 

• the distribution can always be chosen positive, and 

• there are mappings to a second-order differential 
form for all two-body operators. 

From these properties, we show that positive-definite 
Fokker-Planck equations exist for many-body fermionic 
systems, provided that the distribution tails remain suf- 
ficiently bounded. Such Fokker-Planck equations enable 
first-principles stochastic simulation methods, either in 
real time or at finite temperature. As is usual in such 
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methods, care must be taken with sampling errors and 
boundary terms due to the distribution tails. Due to the 
non-uniqueness of the representation, there is a type of 
gauge freedom in the choice of stochastic equation. We 
show how this stochastic gauge freedom, which has been 
successfully used to remove boundary terms in bosonic 
representational, can in principle also be used here. 

Representations for fermionic density operators were 
introduced by Cahill and Glauber 33 using fermionic co- 
herent states^. These provide a means of defining quasi- 
probabilities for fermionic states analogous to the well- 
known bosonic distributions^ 3 ^. However, the result- 
ing quasi-probabilities are functions of non-commuting 
Grassmann variables, and are thus not directly computa- 
tionally accessible. Nevertheless, fermion coherent states 
and Grassmann algebra are useful for deriving analytical 
results in Fermi systems. 

The Gaussian method introduced here overcomes the 
problems inherent in using Grassmann algebra variables. 
The Gaussian expansion utilises an operator basis con- 
structed from pairs of operators, instead of a state- vector 
basis. Because pairs of fermion operators obey com- 
mutation relations rather than anti-commutation rela- 
tions, a natural solution of the anticommutation prob- 
lem is achieved. The resulting phase space thus ex- 
ists on a domain of commuting c-numbers, rather than 
anti-commuting Grassmann variables. Furthermore, the 
phase-space equations obviate the need to evaluate large 
determinants in simulations. This method substantially 
generalizes and extends earlier phase-space techniques 
used in quantum optics to treat electronic transitions 
in atomsSS* 3 ^. It is different to auxiliary field quan- 
tum Monte Carlo methods 3 ^ in condensed matter theory, 
which use Gaussian operators, but involve path integrals 
rather than positive expansions of the density matrix. 

We begin in Sec[n] by defining the Gaussian operator 
basis on which the representation is based, and intro- 
ducing some convenient notations. In Sec lIIII we define 
the Gaussian representation as an expansion in Gaussian 
operators, and then show how the representation estab- 
lishes a novel class of exact Monte-Carlo type methods for 
simulating the real-time dynamics or finite-temperature 
equilibrium of a quantum system. We show how to map 
quantum operator evolution onto a set of stochastic (real 
or complex) differential equations, and give the corre- 
spondences necessary to calculate physical moments. 

Finally, in Sec. [V] we give examples of the applica- 
tion of the method. These are intended to be illustrative 
rather than exhaustive, and further examples and appli- 
cations in greater detail will be given elsewhere. In par- 
ticular, we note that any nonlinear application requires a 
careful analysis of the issues of sampling error and bound- 
ary term behaviour. For simplicity, we focus on the ideal 
Fermi gas, a generic open system master equation and 
the finite temperature Hubbard model, as well as show- 
ing how to apply gauges to modify the drift evolution. 



II. GAUSSIAN OPERATORS 

Before discussing the Gaussian representation, we first 
introduce the fermion operators on which it is based. 
Fermionic Gaussian operators are defined as exponentials 
of quadratic forms in the Fermi annihilation or creation 
operators. This simple definition encompasses a wide 
range of physical applicability. Obviously, it includes 
the well-known thermal density matrices of the free field. 
Since the definition includes quadratic forms involving 
pairs of annihilation or creation operators, it also encom- 
passes the pure-state density matrices that correspond to 
the BCS states used in superconductivity. 

A more subtle issue is that the definition is not re- 
stricted to Hermitian operators. This has the advan- 
tage of leading to completeness properties that are much 
stronger than if the definition were restricted to only Her- 
mitian operators. Some of these issues are discussed else- 
where, in a more formal derivation of the mathematical 
properties of the Gaussian operators^. 



A. Notation 

Before giving mathematical results, we summarize the 
notation that will be used. We can decompose a given 
fermionic system into a set of M orthogonal single- 
particle modes, or orbitals. With each of these modes, 
we associate creation and annihilation operators and 
bj, with anticommutation relations 
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where j, k — 1...M. Thus, b is a column vector of the 
M annihilation operators, and W is a row vector of the 
corresponding creation operators. 

For products of operators, we make use of normal and 
antinormal ordering concepts. Normal ordering, denoted 
by is defined as in the bosonic case, with all anni- 

hilation operators to the right of the creation operators, 
except that each pairwise reordering involved induces a 
sign change, e. g. : bM :— —bjbi ■ The sign changes are 
necessary so that the anticommuting natures of the Fermi 
operators can be accommodated without ambiguity. 

To enable the general Gaussian operator to be written 
in a compact form, we use an extended-vector notation: 



b = 



(2.2) 



is defined as an extended column vector of all 2M oper- 
ators, with an adjoint row vector defined as 



i = (bib 7 



(2.3) 
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Throughout the paper, we print vectors of length M and 
M x M matrices in bold type, and index them where 
necessary with Latin indices: j — 1,...,M. Vectors of 
length 2M we denote with an underline, while 2M x 2M 
matrices are indicated by a double underline. These ex- 
tended vectors and matrices are indexed where necessary 
with Greek indices: /j, = 1, ...,2M. For further examples 
of this notation, se o n i 31 . More general kinds of vectors 
are denoted with an arrow notation: A . 
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B. Definition of the Gaussian operator 

We define a Gaussian operator to be any normally or- 
dered, Gaussian form of annihilation and creation op- 
erators. Like a complex number Gaussian, the opera- 
tor Gaussian is an exponential of a quadratic form, with 
the exponential defined by its series representation. The 
most general Gaussian form is a cumbersome object to 
manipulate, unless products of odd numbers of opera- 
tors are excluded. Fortunately, restricting the set of 
Gaussians to those containing only even products can be 
physically justified on the basis of superselection rules for 
fermions. Because it is constructed from pairs of opera- 
tors, this type of Gaussian operator contains no Grass- 
mann variables. 

With the extended-vector notation, we can write any 
general Gaussian operator A as: 



where '+' denotes a generalised transpose operation de- 
fined by 



(2.8) 



Any Gaussian operator can always be written with a co- 
variance matrix that has the antisymmetry g_ = —g_ + . 

The Gaussian operators are defined here in terms of the 
full 1 + p = 1 + 4M 2 amplitude plus covariance matrix 
components. Alternatively, we can include just the pa- 
rameters A that lead to distinguishable Gaussians, where 



A = (O, n, m, m H 



(2.9) 



A(\) = Q-^ : exp 



-b Zb/2 



(2.4) 



where is an amplitude, TV is a normalizing factor de- 



fined so that Tr 



A(A) 



giving only 1 + p = 1 + M (2M — 1) parameters. We will 
index over the phase-space variables with the notation 
A Q , a — 0, . . .p. For simplicity, we generally deal with 
the full, unconstrained g_ matrices. 

The normalisation A/'contains a Pfaffian whose square 
is equal to the determinant of the matrix. We will show 
that A/" does not appear explicitly in later results. The 
additional variable plays the role of a weighting fac- 
tor in the expansion. This allows us to represent unnor- 
malised density operators like exp(— /3-ff), and to intro- 
duce stochastic gauges that change these relative weight- 
ing factors in order to stabilize trajectories. 



C. Moments 

Just as with classical Gaussian forms, these generalised 
fcrmionic Gaussians are completely characterised by their 
first order moments (to within a weight factor): 



= CI, and £ is a 2M x 2M 



complex matrix. For later identification with physical 
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observables, it proves useful to write £ in the form: 


Tr 




= tlriij , 


£ = (aT 1 - 2f) , (2.5) 


Tr 
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where a is a generalised covariance matrix and J is the 


Tr 







constant matrix is defined as 



I = 



I 
I 



(2.6) 



It is convenient to introduce complex M x M matrices n 
and n = I — n which, as we show later, correspond to nor- 
mal Green's function for particles and holes respectively, 
and two independent antisymmetric complex M x M ma- 
trices m and m+ that correspond to anomalous Green's 
functions. These are related to the covariance matrix a 
by 
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(2.10) 



If the Gaussian operator happens to be a physical den- 
sity matrix, these quantities correspond to the first-order 
correlations or Green's functions. Thus, in many-body 
terminology, n and n are the normal Green's functions 
of particles and holes, respectively, and m and m + arc 
anomalous Green's functions. From this we see that, for 
the subset of Gaussians that are physical density matri- 
ces, we must have that m^m" 1 ", and n^=n. Further- 
more, n and n must be positive semi-definite (because 
< fe) < 1) . 

More generally, the phase-space function 0( A ) corre- 
sponding to the normally ordered operator O is defined 
as a phase-space correspondence, according to: 



(2.7) 



0(A) = (0^_ = Tr 



OA( A) /n 



(2.11) 
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For higher-order moments, a form of Wick's theorem ap- 
plies to any normally ordered product. One simply writes 
down the sum of all distinct factorisations into pairs, with 
a minus sign in front of any product that is an odd per- 
mutation of the original form. The term distinct factor- 
ization means that neither permutation of pair ordering 
nor re-ordering inside a pair is regarded as significant, 
since these do not change the result. Thus an N-th or- 
der correlation (expectation value of a product of 2N 
operators), is the sum of 2A!/(2 Jv iV!) distinct terms, as 
follows, 



: b f j 1 ....b fl2N : 



Here the sum is over all 2Nl/(2 N N\) distinct pair per- 
mutations P(l), P(2N) of 1, 2N, and where (-l) p 
is the parity of the permutation (i.e. the number of 
pair-wise transpositions required to perform the permu- 
tation). 

Thus, for example, the second-order number correla- 
tion moment is: 



(2.13) 



D. Generalised thermal states 

An important subset of the Gaussian operators is the 
set of generalized thermal operators, for which m = 
m + = 0. These include the canonical density matrices 
for free Fermi gases in the case that n, and n are each 
Hermitian and positive definite. More generally, how- 
ever, we do not require n to be Hermitian. In all cases, 
the generalized thermal operators in normally ordered 
Gaussian form can be written most directly in terms of 
the hole population, n = I n: 



A( A) = fidet [n] : exp 



fit 



n" 1 



-2l) T b 



(2.14) 



Of course, there is a symmetry here: in an antinormally 
ordered Gaussian, the role of b' and b is reversed, and 
consequently so is the role of n and n. Our choice of 
normal ordering is in fact arbitrary from a physical point 
of view, and antinormal ordering would also serve our 
purpose equally well, provided all the identities were re- 
defined. 

By comparison, the usual canonical form of the 
fermionic thermal state with a diagonal Hamiltonian 
H = Wuh and a chemical potential p, is an unordered 
form, namely: 



p(r) = exp rh'(p — uj)h jZ . 



(2.15) 



Here Z is the partition function and r = 1/fesT is the 
inverse temperature. In this case, the mean occupation 
numbers are diagonal, and are well-known. They are 
given by the Fermi-Dirac distribution; 



1 + e T ("i-p) 



(2.16) 



However, both Gaussian forms are equivalent. A nor- 
mally ordered thermal Gaussian can always be chosen so 
that n is Hermitian, and hence p(r) = A( A ) if and only 
if il = 1 and mj = rly . 

A rather trivial example is the vacuum state, in which 
n = 0, so that: 



A(l, 0,0,0) 



|0) (0 

: exp — b'b 



(2.17) 



We emphasize that since the Gaussian forms used here 
are not necessarily Hermitian, the generalized thermal 
operators are a much larger set of operators than the 
usual canonical thermal density matrices. 



E. Generalised BCS states 

A second important subset of the Gaussian operators 
is the generalisation of the Bardeen-Cooper-Schrciffcr 
(BCS) states, which are an excellent approximation to 
the ground state of a weakly interacting (BCS) super- 
conductor. The BCS states are the fermionic equiva- 
lent of the squeezed states found in quantum optics, and 
are composed only of correlated fermion pairs. In the 
case of fermions, these are the fundamental pure states 
that carry phase information. In Bose gases, coherent 
states can also carry phase information (as in a laser or 
Bose-Einstein condensate), but the fermionic equivalent 
of these is an unphysical Grassmann coherent state. 

An unnormalized pure BCS state is defined as2& 



|*bcs) = exp [b t gb t /2j |0) , 
so that the corresponding density matrix is: 
Pbcs = *bcs) (*BCs| 



(2.18) 



exp 



b + gb t /2 |0) (0|cxp bg t b/2 



= :cxp btgb^-b^b + bVb^ : .(2.19) 

Apart from being unnormalized, this corresponds directly 
to a Gaussian in our normal form. 

More general non-Hermitian BCS-type states are ob- 
tained on replacing g^ by an independent matrix g + . 
This generalized BCS Gaussian has an extended covari- 
ance matrix of: 



(I + gg+)~ 






(I + g + g) 



-I g 
t+ I 



(2.20) 
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Clearly, from this we can see that the occupation numbers 
and correlations for a generalized BCS state are given by: 



n 
ri 
m 



g+ (I + gg+) g, 
(l + g+g) _1 
(l + gg+) g 
g + (I + ss + ) 1 , 



(2.21) 



which gives the expected result that m + m = fin. 

In summary, the usual BCS states have a density ma- 
trix which is Gaussian, and has g + = g} . These pure 
states exist as a subset of a more general class of BCS-likc 
Gaussian operators. This class also includes operators 
which have g + ^ gt, and are therefore not Hermitian. 
While these operators do not correspond to any physical 
state, a linear combination of them - provided the result 
is Hermitian and positive-definite - can still correspond 
to a possible physical fermionic many-body state. 



III. GAUSSIAN REPRESENTATION 

While the Gaussian operators include a large and inter- 
esting set of physical density operators, there are many 
cases where the existence of interparticle interactions 
leads to more general fermionic states whose correlations 
are of more complex, non-Gaussian forms. In all such 
cases, the overall physical density operator can still be 
expressed as a positive distribution over the Gaussian 
operators. Furthermore, any two-body operator acting 
on a generalised Gaussian can be written as a second- 
order derivative. These important results, proved in 31 , 
means that probabilistic, random sampling methods may 
be used to calculate physical observables, as we show be- 
low. 



A. Definition 

The Gaussian representation for fermion operators is 
defined as an expansion of the density matrix for any 
physical state p{i) as a distribution over the Gaussian 
basis. That is: 



p(t) = J P(X,t)A(X)dX , (3.1) 
where the expansion coefficients are normalised to one: 



P(X,t)dX = 1 . 



(3.2) 



This expansion defines a type of phase-space representa- 
tion of the state: the vector A of Gaussian parameters 
becomes a generalised phase-space coordinate, the func- 
tion P(X ,t) is then a probability distribution function 

over the generalised phase space, and d X = d 2 ( p+1 ^ X is 
the phase-space integration measure. 



B. Moments 

Some basic properties of P(X ,t) follow from those of 
the Gaussian operators. For example, using the normal- 
isation of the Gaussian operators we find that 



Tr[p] = J p(x,t)ridx =n. 



(3.3) 



Thus the normalised distribution P can represent unnor- 
malised density operators by incorporating the normali- 
sation into the mean weight Q. 

More generally, the expectation value of an operator O 
evaluates to 



O 



Tr 



Op /Tr[p] 



P(A,i)Tr 

o(x 



OA 



dx/n 



where the weighted average (. . .) p is defined as: 



0( X 



p(x,t)no(x)dx/n 



(3.4) 



(3.5) 



The phase-space function 0( X ) corresponding to the op- 
erator O is defined as previously, using the generalised 
Wick result of Eq. (|2~T2|) . 

Physical quantities thus correspond to (weighted) mo- 
ments of P. For example, from traces evaluated in 
Sec. Ill Cl we find that the normal and anomalous Green's 
functions correspond to first order moments: 



bib j 



= \Tlij ) p , 
= { m ii) P ■ 



(3.6) 



Number-number correlations correspond to averages of 
products of these moments: 



(: Uifij :) = (nurij 



riijUji + m+m jt ) p 



(3.7) 



where rn = b\bi . 

Similarly, higher-order correlations correspond to 
higher-order moments, the form of which are also de- 
termined by the generalised Wick result of Eq. I|2.12l) . 
We note that the expectation value of any odd product 

of operators must vanish e.g. (bij = 0. Thus the distri- 
bution cannot represent a superposition of states whose 
total number differ by an odd number. Such superpo- 
sition states we exclude from our definition of physical 
state, as they are not generated by evolution under any 
known physical Hamiltonian. The Gaussian distribution 
can, however, represent systems in which particles are 
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coherently added or removed in pairs, leading to nonzero 
anomalous correlations (wiy) p . On the other hand, if the 
total number of particles is conserved or changed only 
via contact with a thermal reservoir, then the anomalous 
correlations will be identically zero and we can represent 
the system via an expansion in only the thermal subset 
of Gaussian operators. 



IV. TIME EVOLUTION 

Here we show how these positive representations of 
density matrices can be put to use. By use of these 
representations, any quantum evolution arising from one 
and two-body interactions can be sampled by classical 
stochastic processes. To see this, note that the time evo- 
lution of a density operator is determined by a master 
equation, of the general form 



d_ 

dt 



p(t) = L[p{t)] 



(4.1) 



where the L is a superoperator that pre- and post- 
multiplies the density operator by combinations of an- 
nihilation and creation operators. 



A. Types of evolution 



We consider three general time-evolution categories: 



Hamiltonian quantum dynamics 

For unitary evolution in real time, the superoperator 
is a commutator with the Hamiltonian: 



L[p\ 



H,p 



(4.2) 



Thermal equilibrium ensemble 

To calculate the canonical thermal equilibrium state at 
temperature T = l/kgr, one can solve an inverse tem- 
perature equation for the unnormalised density operator: 



(4.4) 



the solution of which will generate the unnormalised den- 
sity operator for a grand canonical distribution: p(r) = 
exp[-r(if - fiN)]. 



B. Operator Mappings 

We wish to show how to transform a general operator 
time-evolution equation (Eq. §4.1\i ) into a Fokker-Planck 
equation for the distribution, and hence into a stochastic 
equation. A crucial part of this procedure is to be able to 
transform the operator equations into a differential form. 

The first step is to substitute for p the expansion in 
Eq. EIJ: 



dP{ A , t) 
dt 



A( A )d A 



P(X,t)L 



a(a; 



d\ 



(4.5) 



Second, we use the differential identities derived intfi to 



convert the superoperator L 



A 



into an operator C 



A 



that contains only derivatives of A. Next we integrate by 
parts to obtain, provided that no boundary terms arise, 



dP{ X , t) 
Jt 



A(X)dX 



£ 



P(X,t) A(X)dX 



(4.6) 



where £ is a reordered form of C, with a sign change 
to derivatives of odd order. Finally, we see that this 
equation holds if the distribution function satisfies the 
evolution equation 



Irreversible quantum dynamics 

More generally, for an open quantum system, there will 
be additional terms of Lindblad form22iiiI to describe the 
coupling to the environment: 



H,p] +J2(2d K pd^ - [d K d K ,p\. 



K 



(4.3) 



where the operators Ok depend on the correlations of 
the environment or reservoir, within the Markov approx- 
imation. 



dt 



P(X,t)= £ P(X,t) 



(4.7) 



This procedure for going from the master equation for 
p to the evolution equation for P can be implemented 
using a set of operator mappings, in which we introduce 
antmormal ordering as the opposite of normal ordering, 
and denote it via curly braces: {bjbi} — — b^- . More 
generally, we can define nested orderings, in which the 
outer ordering does not reorder the inner one. For exam- 
ple, {: pbj : hi} — —bib^ : p : , where p is some density 
operator. When ordering products that contain the den- 
sity operator p, we do not change the ordering of p itself; 
the other operators are merely reordered around it. 

Including all possible orderings, we obtain the follow- 
ing mappings: 
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pbb 



— OCT — 



+2a--a 
— —ocr — 



—acr — 



p, 



p. 



;(-) 



— OCT — 



P. 



(4.8) 



Here, er = J — cr, and cn s J = | (er — The notation 

■jr- indicates a differentiation on both left and right sides 

ax G 

with the ordering of matrix multiplication preserved, so 
that: 



(4.9) 



lJ dnik 
d 



d . . 
riiknij 



13 dm k 
d 



%J dni k 
d 



on, 



niknij 
niknij 



P. 
P. 
P, 
P. 



(4.11) 



C. Fokker-Planck equation 

To be able to sample the time evolution of P with 
stochastic phase-space equations, which is the final goal, 
we must have an evolution equation that is in the form 
of a Fokker-Planck equation, containing first and second 
order derivatives: 



dt 



P(X,t) 



a=0 
1 P 
a,b—0 



d d 



D ab {\) 



(4.12) 

P{X,t) , 



For convenience of the reader, these identities are sum- 
marized in a more explicit form using the M x M sub- 
matrices, in the Appendix. We note here that the mixed 
identities involving nested orderings are not independent 
- one can always be obtained from the other. Also, since 
the kernel is analytic, the distinct analytic derivatives 
of the kernel are all interchangeable and lead to equiva- 
lent identities, so that generically if A Q = + iX^, then 
d/dX a — d/dX* = —id/dX v a . Another freedom is that a 
can be replaced by — a + in any of the identities. 

If there are higher than quadratic terms present, the 
differential mappings are applied in sequence. The oper- 
ator set closest to the operator p leads to the innermost 
differential operator acting on P. Thus, for example, 



where a = 0, . . .p is an index that ranges over all the 
variables in the phase space. The matrix D ab must 
be positive-definite when the Fokker-Planck equation is 
written in terms of real variables. Fortunately, the fact 
that the representation kernel A( A ) is analytic in the 

phase-space variables A means that the matrix D ab can 
always be chosen positive-definite after it is divided into 
real and imaginary partseS, through appropriate choices 
of the equivalent analytic forms d/d\ a — d/dX* — 
-id/dXl 

A Monte-Carlo type sampling of Eq. (|4.12() can be re- 
alised by integrating the Ito stochastic equations 

dX a (t) = A a (\)dt + J2 B *bCx)dW b (t), (4.13) 

b 



pbAK'K' 



M 


n 9 


\l' 1/' 




W j 









-Vfi'pO'a.v' 



P 



P (4.10) 



For a system in which the total number is conserved, 
one can use the simpler thermal subset of these corre- 
spondences, i.e. including only those that contain terms 
that remain when all anomalous correlations vanish: 



where dW b {t) are Weiner increments, obeying 
{dW b (t)dW b >{t')) = 5 bM S(t - t')dt, i. e. Gaussian 
white noise. The noise matrix B ab is related to the 
diffusion matrix by D ab = ^ c B ac B bc . This equation is 
directly equivalent to a path-integral in phase-space, so 
that the procedures outlined here can be regarded as a 
route to obtaining a path-integral without Grassmann 
variables. 

Auxiliary field methods— can also be used to obtain 
a non-Grassmann path integral, but these are generally 
much more restrictive. 
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D. Stochastic gauges 

The final phase-space equations are far from being 
unique. This freedom in the final form arises from dif- 
ferent choices that are made at different points in the 
procedure. The choices at some points are constrained 
by the need to generate a genuine Fokker-Planck equa- 
tion with a positive-definite diffusion matrix and vanish- 
ing boundary terms. Other than this, the choices are in 
principle free; they affect the final stochastic behaviour 
without changing observable moments. They are thus a 
stochastic analogue of a gauge choice in field theories, 
and a good choice of stochastic gauge can dramatically 
improve the performance of the simulations^. 

Because the Gaussian basis is analytic, methods pre- 
viously used for the (bosonic) stochastic gauge positive- 
P representation are therefore applicable 3 ^^*^. In the 
fermionic case there are three sources of gauge freedom: 



1. Fermi gauges 

For fermionic systems there is a freedom in the choice 
of operator correspondences, arising from vanishing op- 
erator products; any term involving a square of a fermion 
operator, like a^O, is zero. Terms like this (and products 
of such terms), can be added to the Hamiltonian or Liou- 
ville equation without modifying the density matrix. The 
corresponding additional differential terms may not van- 
ish, hence generating a different but equivalent stochastic 
equation. Such a fermionic stochastic gauge is necessary 
to avoid complex weights in imaginary-time simulations 
of interacting systems, such as the Hubbard mode&. 



2. Diffusion gauges 

Diffusions gauges arise from the fact that the matrix 
square root D ab — J2 C B ac Bbc has multiple solutions, es- 
pecially if one notes that there is no restriction on the sec- 
ond dimension of B ab . This changes the stochastic noise 
term and can lead to a reduction in sampling erron4^. 



3. Drift gauge 

Drift gauges are obtained by trading off trajectory 
weight against trajectory direction. The possibility for 
drift gauges arises from the weight f2 in the density- 
operator expansion. The first of the correspondences in 
Eq. (|4.8|) can be used to convert drift terms for the phase- 
space variables into diffusion terms for the wei^hlii^. As 
a result, one can add an arbitrary gauge g a ( A ), of the 
same dimension as the noise vector. Assuming E>o b = 0, 
and using Einstein summation conventions, this leads to: 



dfi(t) = A dt + Qg b dW b (t), (4.14) 
dX a (t) = A a dt + B ab [dW b (t)-g b dt]. 

Previous work^Siii has shown that drift gauges can re- 
move boundary terms in bosonic positive-P representa- 
tion by stabilizing deterministic trajectories. 

V. EXAMPLES 

The virtue of phase-space representation is that while 
Hilbert space dimension grows exponentially with the 
number of modes M, the phase-space dimension only 
grows quadratically. Thus, for example, a problem in- 
volving M = 1000 fermion modes has a Hilbert space di- 
mension of D — 2 1000 = 10 103 dimensions. This is larger 
than the number of particles in the observable universe 
(which is perhaps 10 85 by current astrophysical reckon- 
ing). By contrast, the fermion phase-space dimension is 
4 x I0 6 . While large, this is not astronomical. 

Hamiltonians and general time-evolution equations 
that are only quadratic in the Fermi ladder operators, 
i. e. constructed from one-body operators, will map to 
a Fokker-Planck equation that contains only first order 
derivatives. The evolving quantum state can thus be 
sampled by a single, deterministic trajectory. More gen- 
erally, quartic terms and cubic terms (if Bose operators 
are included) can also be handled, and these result in 
stochastic equations or their equivalent path integrals. 

Examples of how some typical Fermi problems are 
mapped into phase-space equations are given as follows. 

A. Free gas 

As an example of quadratic evolution, consider the 
thermal equilibrium calculation for a gas of noninteract- 
ing particles. The governing Hamiltonian (including the 
chemical potential) is always diagonalizable, and can be 
written as: 

H = b t u:b, (5.1) 

where u, 3 = SijLOj&Te the single-particle energies. The 
grand canonical distribution at temperature T = I/fc^r 
is found from the equation 

^-p = -i(btu,bp + pbtu;b) . (5.2) 

Now this master equation can be mapped to an equiv- 
alent equation for the the distribution P by use of the 
thermal correspondences in Eq. QJ. However, because 
the solution is an unnormalised density operator, there 
will be zeroth-order terms in the equation. We can con- 
vert such terms to first order by applying the weight (f2) 
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identity in Eq. (|4.8|) . thus obtaining the Fokker-Planck 
equation 



dP 



k 



d d ' 

an k oil 



n k P . (5.3) 



This Fokker-Planck equation with first-order derivatives 
corresponds to deterministic characteristic equations: 



17 = - uj k ttn k , 



-u k n k (l-n k ) 



(5.4) 
(5.5) 



Integrating the deterministic equation for the mode oc- 
cupation n k leads to the usual Fermi-Dirac distribution: 



1 



n k 



(5.6) 



From integration of the weight equation, one finds that 
that normalisation of the density operator is 



Tr [p u ] = n(r) = n n k e-^ n " T , 



(5.7) 



i. e. the weight decays exponentially at a rate given by 
the total energy. 



B. General quadratic evolution 

More generally one can have a quadratic Liouville op- 
erator in situations involving non-thermal terms like bibj. 
This can occur for, example, when fermion pairs are 
generated from molecule or exciton dissociation. These 
are even associated with certain spin-chain problems 43 , 
where the Jordan- Wigner theorem is used to transform 
spins to fermion operators. Other quadratic Liouville op- 
erators are commonly found in cases involving coupling 
to reservoirs^. 

The generic phase-space equations for a general Fermi 
system with a quadratic Liouville operator can be easily 
obtained, for evolution both through time and through 
inverse temperature. The most general master equation 
that covers both kinds of evolution can be written 

C„ : +C; u {:pt fl :bi}) , (5.8) 

where the elements of 2M x 2M matrices A, B and C 
are determined by the coefficients of the Hamiltonian or 
master equation. By applying the mappings of Eq. I|4.8|l . 
we find the evolution of the covariance matrix to be: 



eh 



a = aAa + aBa + aCa + aC^a 



(5.9) 



This equation simply corresponds to the characteristic or 
drift equations given by the vector A in the Ito stochastic 
equation l|4.13[l , and in these cases there is no diffusion or 
stochastic term. Unlike a conventional path integral, we 
see that a quadratic Hamiltonian or Liouville equation 
simply results in a noise-free, deterministic trajectory on 
phase space. For deterministic evolution such as this, the 
weight SI does not affect physical observables, so we do 
not consider it here. 

In the examples that follow, we assume for simplicity 
(but without loss of generality) that the constant matri- 
ces have been chosen with hermitian anti-symmetry such 
that: 



A = -A+ 
jg = -£+ 
g} = -g+ . 

Temperature evolution 



(5.10) 



For temperature evolution, the structure of the master 
equation (Eq. |4.4|l) is such that A — B_ and C = Cj , 
giving the simpler result: 



1 



dr= 2 

where we have introduced: 



(I -2a) T (I -2a) + a 



(5.11) 



T 



g-g 
1 



(5.12) 



For the case of a number conserving Hamiltonian H 
b^wb, we find that B_ — and 



-uj t 
u 



The phase-space equations then reduce to 



1 



(nam + nun) . 



(5.13) 



(5.14) 



dr 2 
which reproduces the free gas example above. 

2. Dynamical evolution 

For time evolution, with possible coupling to the en- 
vironment, there is a different symmetry to the master 
equation (Eq. l|OJ|) that means that A + B-C-C} = 0. 
A formal solution to the phase-space equations can now 
be explicitly written down: 

gSt) = exp (g(0) - g°°) exp (-gt) + g°°, 

(5.15) 
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where U_ = {B — Q)I_ and where a°° satisfies: 

LgL = V} g°° +g°°IL- (5.16) 

To illustrate the physical meaning of these matrices, 
we consider the simplest model of a small quantum dot 
coupled to a zero-temperature reservoir: 



iujb^bp + iu'fh 



h^b + 7 (bpV - -6% - \$oj ■ 

(5.17) 



In terms of the general form, this corresponds to A = 0, 
3. = iL and 



—ILO 







(5.18) 



The general solution then reduces to 



s(t) 



,-iu>--y/2 



Juj — y/2 







Ju+j/2 



ji(o)-I) 
-I, 



(5.19) 



which implies that the density decays as n(t) 



-it 



n(0), 



as expected. 

The solution to a multimode quantum dot model also 
follows from Eq. I|5.15|l . The relevant master equation is 

p = -iujibfijp + itJjipTffij 



hpb 



' 2 

for which the evolution matrix is 

D -iu T +-y T /2 



l $biP- \pb^ , (5.20) 



U = 









Ju+~y/2 



(5.21) 



Physically, this corresponds, as expected, to damped os- 
cillatory behavior (taking 7 to be positive definite) in the 
moments: 



n = e 



J— Y/2 



n(0)e 



-iu> — y/2 



'-l T /2 



m(0)e 



—iw — y/2 



(5.22) 



Here, of course, there are no electron-electron interac- 
tions included. However, such interactions can be dealt 
with via a stochastic sampling methods, as we show in 
the next section. 



C. Interacting gas 

1. Two-body interactions 

For systems of particles with two-body interactions, 
the Gaussian representation gives nonlinear, stochastic 



phase-space equations, which must be solved numerically. 
Consider a two-body interaction of the form: 



H 2 = 2J UijUifij , 



where ru 



a a.j 



For a number-conserving system, 



we can use correspondences of Eq. (I4.11|) to generate a 
Fokker-Planck equation for the grand canonical evolu- 
tion. The diffusion matrix D UiV in this equation is 



D 



ij,kl 



pq \WipTlpjTlkqTlqi 



TlipflpjUj^qTlqiy . 



(5.23) 



Suppose that the interaction matrix U pq is negative- 
definite, such that we can write it as a sum of negative 
squares: U pq = — ^ a b p ^ a b q ^ a . Then the diffusion matrix 
is positive definite, as it can be written in the form: 



D 



ij.kl 



E 



) B w 



y. ( 5 - 24 ) 



where the noise matrices are: 



B 



(2) 



^ ^ bp.aTlipflpj 
P 



(5.25) 



Thus for an interaction of this type, the noise terms in 
the final stochastic equations will be real. The form of 
noise terms for a more general interaction is considered 



,46 



2. Hubbard model 

As an example, we apply the representation to 
the Hubbard model, which is the simplest nontriv- 
ial model for strongly interacting fermions on a lat- 
tice. It is an important system in condensed matter 
physics, with relevance to the theory of high-temperature 
superconductors^ 6 , and in ultracold atomic physics. The 
full phase-diagram in two dimensions is not fully under- 
stood as yet. Due to developments in atomic lattices, 
this model is directly experimentally accessibl o 6 i 44 . 

The Hamiltonian for the model is2&: 

fr(ni,n_i) = E " ' r E " 1 "<••<• ; - 



where n,- 



5-t 



(5.26) 

;{n CT } - - . The index a denotes spin 



HJ,<7 — "-i lCT "-J," L"»Jtj- 

(±1), the indices i,j label lattice location. Here ty = t 
if the i,j correspond to nearest neighbour sites, = p 
if i = j and is otherwise 0. The chemical potential p is 
included to control the total particle number. 
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Because the Hubbard model conserves total number 
and spin, one can map this problem to a reduced phase 
space of A = ($1, ny^n^-i). Thus the simpler mappings 
of Eq. (|4.11l) can be used for each spin component. The 
one-body terms generate drift terms only, and can be 
dealt with as above. The two-body terms generate both 
drift and diffusion terms. Applying the mappings directly 
to the Hubbard model as written above, we obtain the 
diffusion matrix 

Dija^klfj' U3<j,— a' ^ {TljpgTlpjcrTlkpu' Tlpl a ' 

P 

^Hipcr^pja^kpcr'^pla'} j (5.2T) 

which, because it has zeros on the diagonal, cannot be 
put into a positive definite form with real variables. 

However, using the anticommuting properties of the 
Fermi operators, we can rewrite the interaction term in 
the Hubbard Hamiltonian as 

Hi = -- ^- : fei _ S^jj-i) 2 ■ 

3 

— ^ Uifj^ja 1 • (5.28) 
3 

where S = U/\U\ = ±1. Now in this form, the interaction 
matrix is negative definite: 

Uiaja' ~ Sij ^cT,a' S5a. — a'^j 

k 

where s = (S + l)/2, so that s — for the attractive case 
and s — 1 for the repulsive case. 

From Eq. (|5.24(l the diffusion matrix is positive defi- 
nite, with corresponding noise matrices: 

B§l,a = V\U\j2v S n la n aj , 

B§l, a = VW\]2o s n la n ao . (5.30) 

With this choice of noise terms, the final phase-space 
equations are, in Ito form, 

^ = ijn.rWn. + n^W}, (5.31) 

where we have introduced the stochastic propagation ma- 
trix: 

= Ui ~ -\, {' » ., a + <r"$ r) } • (5-32) 

The real Gaussian noise Q r ' (r) is defined by the correla- 
tions 

(£j r) (r)#V)) = 2\U\5(r-r')6 jf S rr/ . 

Because the diffusion can be realised in terms of real 
noise, the phase-space equations will not be driven off 



the real manifold. This has an important implication 
for the weight O, which enters the problem because the 
solution will be an unnormalised density operator. The 
weights for each trajectory evolve as physically expected 
for energy-weighted averages, with weights depending ex- 
ponentially on the inverse temperature r and the effective 
trajectory Hamiltonian H: 

^=-JlH(n 1 ,n_ 1 ) . 

(XT 

Because the equations for the phase-space variables rii^ a 
are all real, the weights will all remain positive, thereby 
eliminating the traditional manifestation of the sign prob- 
lem. 

This method can calculate any correlation function, 
at any temperature, to the precision allowed by the sam- 
pling error and subject to there being no boundary terms 
in Eq. (|4.6|) . Preliminary simulations in one^, twos* 
and three dimensions showed that sampling error is well- 
controlled, even for very low temperatures. However, 
more extensive simulations of the 2D Hubbard model 
have shown that, at half filling, certain correlation func- 
tions do not appear to converge to the correct zero- 
temperature results at these very low temperatures^. 
Because the Gaussian basis does not possess many of the 
symmetries of the Hubbard model, they must be restored 
in the distribution over Gaussian basis elements. For fi- 
nite sampling, this restoration may be incomplete, giving 
the departure from exact results at low temperatures. It 
has been shown that the correct results can be obtained 
by applying a projection onto a symmetric subspace^. 
There may also be systematic errors if boundary terms 
are present. Both of these possibilities imply that fur- 
ther optimization via stochastic gauge choices may be 
required to keep the low-temperature distributions com- 
pact, free from tails and from features that would lead to 
biasing. 

3. Drift gauges 

For the Hubbard model, we can modify the drift part 
according to Eq. 14.1411 by adding a term G a to the 

(r) 

stochastic propagation matrices . Because of the di- 
agonal nature of the noise terms, the added term will 
also be diagonal: Gy 1<T = <%Crj )CT . The additional diffu- 
sion term in the weight equation is then 

(f ) = wi^' G " ( ' K (5 ' 33) 

x ' 9 1 jra 

The choice of gauge term Go- is guided on the one hand 
by the need to ensure the phase-space distribution re- 
mains bounded and on the other by the requirement of 
introducing only the minimum amount of diffusion into 
the weight. The function should thus act only when nec- 
essary to control large trajectories and should be zero 
otherwise. 



12 



The diagonal form of gauge term possible is able to 
remove instabilities in the Hubbard equations that are 
directly due to the interaction term U. However, insta- 
bilities may still arise from the coupling terms tij, even 
though they are of lower order. Thus it may be neces- 
sary to introduce weaker, off-diagonal gauge terms. This 
in turn requires additional, off-diagonal noise terms in 
the propagation matrix. Such noises can be introduced 
by use of additional Fermi gauges. For example, the van- 
ishing term^I 

(5.34) 

where Vij jCT are positive numbers, gives the additional 
stochastic contribution to the propagation matrix: 



T, 



(r) 



(r) 



r£Uc8i(T) 



(5.35) 



where the new noises Qj a (r) have the correlations 



4 r) (r)4 r P(r')) = 4^,^(T-r')fe^M5.36) 



We can now introduce arbitrary off-diagonal gauge terms 
Gi.^ a into the propagation matrix, with the correspond- 
ing diffusion term in the weight equation 



^E G y>4 r / 4 ^>- ( 5 - 37 ) 



Again there is a trade-off between gauge strength and 
additional diffusion. But there is also a freedom (in the 
choice of Vij^a) as to whether the noise appears in the 
weight equation or in the propagation matrix. 

With such a combination of Fermi and drift gauges, it 
is possible to introduce terms to stabilise the drift evo- 
lution of any of the phase-space variables riy. CT , and so 
maintain a bounded phase-space distribution. 



of application. We have given examples of the use of 
fermionic differential identities to transform multi-mode 
master equations into deterministic phase-space equa- 
tions, although more general interactions typically lead to 
stochastic equations. These equations have exponentially 
less complexity than the full Hilbert space equations, are 
generally simpler to solve than path integrals, and never 
involve either Grassmann variables or determinants. 

The application to the Hubbard model demonstrates 
the immediate utility of the Gaussian method to solving 
long-standing problems in many-body quantum physics, 
provided suitable gauges can be found to ensure that 
boundary terms to not arise. Rapid experimental ad- 
vances in the area of ultra-cold fermionic atoms^ mean 
that direct and quantitative tests of precise theoretical 
predictions should be feasible in the near future. Demon- 
stration of a quantum degenerate Fermi gas in a lattice 
has already taken placed. 

The general technique established here potentially also 
has broad applicability in many other areas of quantum 
many-body theory and quantum field theory. 
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VI. CONCLUSION 



In summary, we have introduced a phase-space repre- 
sentation for many-body fermionic states, enabling new 
types of first-principles calculations and simulations of 
highly correlated systems. Many-body systems with 
one- and two-body interactions can be solved by use of 
stochastic sampling methods, since they can be trans- 
formed into a second-order Fokker-Planck equation, pro- 
vided a suitable stochastic gauge is chosen to ensure that 
the distribution remains sufficiently bounded. 

These techniques are potentially applicable to a wide 
range of fermionic problems, including both real-time 
and finite temperature calculations. Generalized mas- 
ter equations for non-equilibrium fermionic open sys- 
tems coupled to reservoirs are a particularly suitable type 
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APPENDIX 

It is sometimes more convenient to work with explicit 
n, m, m + submatrices rather than the total covariance. 



In fully indexed notation, using the M x M submatrices, 
the Fermi operator correspondences (Eq. (^2J))become: 



bibjp 
pbibj 
bjpbi 



pW\ 



d d d 

{nijh lk + m+m jk } - g— {mijh ik + ham jk } + -g—f {ni 3 m+ + m+n kj } 



d d d 

{nijtiik + m+m ]k } + ■ g -— {m lo n ik + nu<m jk } - — - {hijmf k + m+h kj } 



dni k 
d 



dm 



. h {nijfi ik - m+mjk} + g^ {m l3 h lk + h u m jk } + {n y m+ + m+h kj } 



13 dni k 
d 



dnik 

d 
dnik 

d 



{numjk - ni 3 m lk } 
{nijiTiik - humjk} 
{humjk + ni 3 m lk } 







dmi k 

d 
dmi k 

d 
dmik 

d 



{mumjk - mijmik} - - — t- {ni 3 n k i - nun k j} 



dm tk 

d 



{mumjk - mijmik} ~ - — r {hijhki - huh k j} 



dm 
d 



ik 



{mumjk - mijm ik } - ——r {h H nkj - h^nki} 



'Ik 



{m+h ik - m+hjk} - [h^k - h u hjk} - g-j- {m+m+ - m+m+ } 



d 



[m^rijk - mfjiiikj 



m ij - 1 - //',••/',/. r - -7— {njiTHk -nurijk} 



d 
dmik 
d 



dm ik 



{ m U m tk ~ m ki m tk} 



m ^ + dn~^ [ m tj n ih + mfihjk} - g— {nuhjk - h jinik } - — p {m+m+ - m+m+ } 



dm ik 



P, 
P, 
P, 
P, 
P, 
P- 
P, 
P, 
P, 

P, (1) 



where we have used the Einstein summation conven- ference between these equations and the full covariance 

tion for repeated indices. Furthermore, we have explicitly equations is due to the fact that these equations corre- 

written out the extra derivative terms that arise from the spond to a covariance which is constrained to satisfy the 

antisymmetry of m and m + , such that the summation Hermitian anti-symmetry condition, 
of these terms is only for k > I. The factor of two dif- 
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